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ABSTRACT 


This volume is the documentation of the theoretical development of the 
forced steady state analysis of the structural dynamic response of a turbine 
engine having nonlinear connecting elements. Based on modal synthesis, and the 
principle of harmonic balance, the governing relations are the compatibility of 
displacements at the nonlinear connecting elements. There are four 
displacement compatibility equations at each nonlinear connection, which are 
sol ved by iteration for the principal harmonic of the excitation frequency. 

The resulting computer program, TETRA 2, combines the original TETRA 
transient analysis (with flexible bladed disk) with the steady state 
capability. A more versatile nonlinear rub or bearing element which contains a 
hardening (or softening) spring, with or without deadband, is also 
incorporated. 
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1.0 SUMMARY 


The NASA-Lewis Blade Loss program which was originally written for the 
calculation of transient response of turbine engine structures, has been 
extended to predict steady state vibratory response. The original name of 
TETRA (Turbine Engine Transient Response Analysis) has been kept, but slightly 
modified to TETRA 2. This new program can be used to calculate transient and 
forced steady state responses. The 'transient' capability is based on the 
latest version of TETRA and includes the subsequent additions of: 

1) Flexible Bladed Disk module 

2) Squeeze Film Bearing Module 

The 'steady state' capability includes the original modules with the 
exception of the squeeze film module. However, the original rub element has 
been generalized to allow either dead-band or no dead-band, and the rub spring 
characteristics is a hardening spring and/or a linear spring. 

Inputs for transient and steady state are the same except for the obvious 
differences: 


Transient requires 'time-sweep' inputs 
Steady state requires ' frequency- sweep' inputs. 

The basic theoretical approach for the steady state capability requires the 
formation of a global matrix equation in terms of the generalized coordinates 
and nonlinear physical forces. Solution is by harmonic balance and iteration 
of physical displacements at the nonlinear connecting elements. This solution 
yields only the first harmonic of the forcing frequency. 
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2.0 INTRODUCTION 


In the NASA-Lewis sponsored Turbine Engine Transient Analysis program 
(TETRA) , a computational tool was developed to predict the transient dynamic 
response of engineering structures to suddenly applied loads, such as from 
the loss of a blade (1)*. The capability of this program was further enhance 
by the addition of two modules; 1) Flexible Bladed Disk (2), and 2) Squeeze 
Film Bearing (3). The latter was added by Case Western Reserve University 
under NASA-Lewis sponsorship. 

The fundamental technical approach is the method of model synthesis (4), 
wherein the dynamic response of a complex structure is constructed in terms 
of the natural modes of its principal structural components. The possible 
breakdown of a turbine engine into its main components is shown in Figure 
1-1. The equations of motion in the modal generalized coordinates are solved 
numerically by central difference integration in the time domain. This 
solution has the flexibility to accommodate nonlinearities , such as tip 
rubs, squeeze films or other nonlinear bearings or connecting elements. Also 
the gyroscopic coupling between motions in the vertical and horizontal planes 
of rotating structures is considered for rigid as well as flexible bladed 
disks. Applications of the TETRA are found in References 1, 2 and 5. 

The transient response of a structure is a history of the motion and 
loads which initially vary non-uniformly in time until a steady state condi- 
tion is reached. Where damping is low and modal frequencies high, the time 
steps required to reach steady state can be considerable, for this reason a 
more direct method to calculated steady state response was undertaken. 

Steady state capability allows the calculations of forced response 
amplitudes as function of excitation frequency so that engine response from a 
sinusoidal input, such as unbalance, can be obtained over the engine 
operating speed range. For purely linear systems, the methodology is well 
established. 

However, in the presence of nonlinearities, obtaining the steady state 
solution is neither simple nor straight-forward. This is especially true 
with large systems of equations with strong nonlinearities. To date, there 
is no mathematical method to solve the general nonlinear differential 
equations, (6), (7) and (8). 


Numbers in parentheses indicate references. 
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Figure 1-1. Structural Breakdown for TETRA Analysis 


To produce a more pragmatic solution to this complex and important problem, 
the solution is limited to the first harmonic of the forcing frequency. 

Solution by iteration of the compatibility conditions at only the nonlinear 
connecting elements results in a dramatic reduction in the number of equations 
to be solved. The method of harmonic balance, due to Kryloff and Bogoliuboff 
(9), was used to transform the nonlinear differential equations to a system of 
nonlinear connections. 

This methodology was implemented in a computer code built on the original 
TETRA program. The new program, TETRA 2, has the capability of the original 
transient analysis as well as the steady state solution. To make the TETRA 2 
user friendly, the inputs of the original TETRA has been kept unchanged as much 
as possible. The only major addition lies in the description of the 
'time-sweep' of the transient analysis and the 'frequency-sweep' of the steady 
state. The basic description of the structural subsystems and their modal 
data, the concatenation or assemblage of the subsystems and the connecting 
elements are essentially the same. 

This final report is the documentation of the development of TETRA 2, in 
particular, the steady state capability. This report consists of two volumes: 

1) theory, and 2) user's manual. The latter includes the entire input 
sheets for the TETRA 2 code: transient and steady state, as well as sample 

cases, and comparisons of results made for the original two subsystem 
demonstration cases. 

With the NASA sponsored TETRA 2 computer code, industry is presented a 
comprehensive turbine engine rotor dynamics computer code that can be used 
for the calculation of both transient and steady state responses. The nonlinear 
capability of the program greatly enhances and broadens its application to more 
realistic analysis of real engines. 

The authors wish to acknowledge the technical help provided by their 
colleagues. Dr. J. K. Casey contributed his mathematical expertise in the 
computational strategy employed in the program, while M. J. Stallone provided 
overall technical guidance, especially insisting in making TETRA 2 inputs and 
outputs as similar as possible to the original transient version. Also, we 
recognize A. Storace for contribution in the overall program and R. Holt for 
writing the model program that checked out TETRA 2. Finally, to Jeanette Sturgill 
for a very meticulous and neat typing of the working equations; many thanks. 

To our colleagues at NASA-Lewis, thanks to Gerry Brown who managed the 
original TETRA and Bob Kielb who succeeded him in TETRA 2 and Chuck Lawrence, 
NASA-Lewis' monitor for the program. We would like to recognize the thorough 
critical evaluations that the first project monitor of the NASA Blade Loss 
Program, Ming Tang, has made to TETRA and for the subsequent definition of the 
early stages of TETRA 2. He has since passed away last year. His thoroughness 
in reviewing the technical documentations of the Blade Loss Program as well as 
the very cogent constructive criticisms he provided, have proven invaluable. 
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3.0 PART I: ANALYTICAL DEVELOPMENT 


3.1 Technical Approach to Steady State Response Analysis 

The steady state response of systems with nonlinearities is a relatively 
undeveloped field unlike purely linear systems. Nonlinear differential 
equations have many possible solutions, each being highly sensitive to initial 
conditions, external forces and system parameters. For instance, in the case 
of Buffing’s or van der Pol’s equation, grossly different solutions are 
produced by changes in initial conditions or force amplitudes. One initial 
condition may result in periodic solutions; another may result in aperiodic 
solutions; still others may yield jumps or bi-stable solutions, limit cycles, 
and subharmonic and superharmonic oscillations. 

When one considers a complex turbine engine with many degrees of freedom 
and several nonlinearities, the problem of finding general steady state 
solutions is considerable. In most cases, only approximate solutions with 
narrow constraints are practicable. In the present work, a pragmatic 
philosophy is employed which will limit solutions to the first harmonic of the 
excitation frequency and simple harmonic oscillations with constant amplitude. 
This is justified from most engine experience, where engine response is 
essentially in the excitation frequency. 

The use of modal synthesis on a large system also limits nonlinearities to 
inter- subsystem connecting elements, such that the deflection of a structural 
subsystem can be represented by a superposition of its normal modes. The 
nonlinear connecting forces are written as function of the relative physical 
displacements between the joined components and are treated as quasi -external 
forces. By multiplying these connecting physical forces by the appropriate 
modal displacements, the global equations in the generalized modal coordinates 
of the complete assembled structure are obtained. 

The left-hand side of the global matrix differential equation contain the 
linear part which are proportional to the generalized coordinates. On the 
right-hand side are the generalized forces which are functions of time in the 
case of external forces, and the nonlinear connecting forces which are 
explicit functions of the relative physical displacements. 

Employing the principles of 'harmonic averaging' or 'harmonic balance' due 
to Kryloff and Bogoliuboff (9,10) and slowly varying functions (10), the 
differential equations are transformed to nonlinear algebraic equations. After 
solving explicitly (by inversion) for the generalized coordinates, the relative 
physical displacements at the nonlinear connections are obtained by modal 
superposition. These last equations relate the vector of relative physical 
displacements at the nonlinear connections on the left-hand side to the 
external and nonlinear connecting physical forces (in terms of relative 
physical displacements) on the right-hand side. Thus, the equations to be 
solved are reduced from twice the number of generalized coordinates to just 
four times the number of nonlinear connecting elements, which is a considerable 
reduction. 
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Numerical solution by iteration of these compatibility relations is made 
with existing computer subroutines such as those in the IMSL library. 

The following sub-sections delineate the technical approach described here, 
in more technical details. 


3.1.1 Generalized Global Equations of Structural Systems With Nonlinear 
Connecting Elements 

Following the method of modal synthesis, the global equations of 
motion of a system of inter-connected elastic bodies were derived in 
Reference 1,2, and 4. Those obtained in (1) and (2) were solved in 
the transient TETRA and are the same ones solved in TETRA 2. To 
illustrate the basic concept of the method employed in the steady 
state analysis, the denotation of each structural component and 
plane(s) (horizontal and vertical planes) will be omitted for 
simplicity. A very detailed development of the working equations 
actually programmed, is presented in a later section. 

Conceptually, the equations of motion of the fully coupled and 
assembled system are written in tensor (or matrix) notation. The 
dependent variables, q^(t) are the generalized modal coordinates of 
the component structures. Since each component structure’s normal 
modes are obtained with free-free boundary conditions, these satisfy 
orthogonality only within the individual structures, such that 
coupling forces between substructures and their normal modes exist. 

In short, the matrices of the coefficients of the generalized modal 
coordinates are not diagonal, but in general, full matrices. 

The generalized nonlinear connecting forces are functions of the 
relative physical displacements between the joined structures at the 
connected points. Though these displacements belong to a subset of 
the physical displacements of the substructures, these ’connection’ 
displacements are not expressed in terms of normal modes. Thus, the 
nonlinear connecting forces are always in terms of the relative 
physical displacements. 

From (1) and (2), the global equations of motion in the 
generalized coordinates are given in tensor form as: 


Where: 


modal mass matrix (non-diagonal with the flexible bladed disk 
module) 

“ij = modal natural frequencies, diagonal matrix 

^ij - modal critical damping ratio 

Gij= gyi^o coupling - skew symmetric matrix 

n » axial location of linear connecting element 

» generalized linear connecting spring matrix, summed over 
the ’n’ linear connections 

Cj j'^=general ized linear damping matrix associated with 

F^(t)- generalized external force - sinusoidal 

= modal nonlinear connecting forces at m locations, 
i Summed over m locations. 

Ax”** Relative physical displacement at the nonlinear connecting 
element with axial location m. 

The generalized nonlinear connecting forces are given as a 
hardening spring with viscous damping, and written in terms of the 
physical displacements. For example, one nonlinear spring force at 
point ’m’ is given as a physical force: 

F"* = 


Where: K"** the linear part of the spring 

m 

^ - the nonlinear factor 


X"** the relative physical displacement at location 'm’ 
viscous damping coefficient of the connecting element. 


The generalized nonlinear connecting force on a subsystem at 
point ’m’ is simply the product of the physical force F*^ times the 
subsystem’s modal displacement at point ’m’ in the i^" mode; thus: 


(h”*F'" = - cj)'" 


m.2 


1 + p'"(AX ”) 


Where: 


AX'" + o'" AX" 




= I 


:th 


mode displacement at point 'm' 
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When there is more than one nonlinear connection, the physical 
connecting forces are multiplied by the modal displacements at their 
locations and summed over those locations - for each mode. The i^^ 
generalized nonlinear connecting force is therefore: 


H = -Y <D'"F'"= - y o'" 

I C i 




Ax" 


- 2 <t>'"C'"AX" 


The global generalized equations may therefore be rewritten: 


'V. = M..q.+ 

I ipj 




n ^ 

- Z 1 AX'" - X <i)'"C'"AX'" = 0 

m 'h 

3.1.2 Transformation bv Harmonic Balance and the Compatibility Relations 


The global generalized differential equations are nonlinear in the 
physical displacements at the nonlinear connecting elements. To obtain 
forced solutions synchronous with the excitation frequency, an 
assumption of slowly varying function is made. This assumes that the 
principal motion is harmonic with constant amplitude, so that the 
ensuing generalized - and physical - displacements, may be written in 
the form: 

Q U) = a .cos(jit + b sinoit 
X'"(rt = A'"cos«i + 

with i. = 6. = A'" = S'" = o. = 6. = A'" = S'" = 0 

(I t ( 

The next step is to employ the principle of harmonic balance 
(9,10,11) by means of which the "N" differential equations are 
transformed to "2N" nonlinear algebraic equations.* 

The transformation is made as follows: Substitute the harmonic form 

of the displacements into the N global generalized differential 
equations; multiply each equation by 'coswt dt’ and integrate from 0 
to 2n/o) ; repeat the last step but multiply with 'Siruat dt’ and 

integrate as before. 


Recalling that each global generalized equation is identified as: 
'Vi = 0, the transformation becomes: 



2n/(o 

0 


'V, 


a 6 A'"B'"&)f... 

V I 



dt=0 


Where: 'Vi - cosine transformation 

'V- = sine transformation 


*2N, when motions are in only one plane; 4N when motions are in 2 planes. 
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The following integral transformations, which are inferred from the 
terms of the global generalized differential equations, are required: 



0) 

n 


■ 2n/o) 


[X 


' t=o 


mi3 / cos(j)t 


SlfKjit 




dt = 


3 m 
- A 

4 

3^m 

4 






Where: 




a.cosbit + bsinuit 

I I 


X"* = A"*cos(ot + B"^sinu)t 


With these transformations, the differential equations become 
algebraic equations. Thus, 


AA.. 

y 

BA.. 

y 



ij: 1,2 ... N 


Where: AA.jj, AB,jj, BA,-j, BB.,-,- are the coefficients of 

“^the generalized coordinates 

- cosine transformation of the external general force 

■ Sine transformation of the external general force 

H": - Cosine transformation of the nonlinear connecting force, a 
function of the A*", B™ components at all nonlinear connecting 
elements. This may be a polynomial in A"’, B"’ of the form: 

(A‘)^(B‘)^(A‘)(A')^«c. 

Where m: i,2,...Af 

■ = Same as except it is for the Sine transformation. 

The previous transformed matrix equation may be written more 
simply with a change in the range of the index i: 


li 


AA"*, AB"*... 


Where 


E = F . + H 

ij j ‘ I 

iJ = 1,2,. ..2N 
ct. = a., i=\,2,...N 

«^, = h., ,N+2,...2N 


E . = 
y 


AA 


Inverting the Eij matrix, one obtains an explicit expression for the 
generalized coordinates: 


C. = y ( AA"*,AB"*... 

* iJ J y 7 [ ■" \ / J denotes 


index j 
summation 


The first term on the right-hand side is the vector of generalized 
displacements due to the external forces alone, while the second are due 
to the unknown physical displacements at nonlinear elements. 


11 


Thus there are 2N equations with 2N + 2M unknowns; the latter being 
the cosine (A"*) and Sine (B"*) components of the relative physical 
displacements atAf nonlinear connections. These are the ultimate 
unknowns. To solve for the A"* and B'^’s, displacement compatibility 
relations at the nonlinear connections are formed. 

Recall that: = 


From the assumed harmonic form of the displacements, it follows that 
at point m: 

A" = fr»i = i = I 

i=l i=l 


2N 


B- = = X ♦r*, = I 


i = l 


i=AT+l 


where the repeated index i indicates summation over the modes. 

The components of the physical displacements at any of the nonlinear 
elements are therefore: 


NOTE: 



N 

n 

r ^ 

r 

A'" 

= I = 

■■ i-f," 


\ II 1 


i = l 

j=i 

‘y-i 

[ -J J 


2N 

2n 


r 2N 

B"' : 

= z 4>r«, 

= Z 

1 

Y 


i=iV+l 

i=iV+l ' 

j = N + l 


N 

2N 




- <}>. 





i=l ‘ 

i = N + i 






Forming similar physical displacements at both ends of a nonlinear 
connecting element, one obtains their difference which is the relative 
physical displacement. Thus, for the cosine component at point m; 


N 


®a'”-©a'" = aa' 




j = l 


N 

1 

J = 1 




+ E~^H 


‘J J\ — 


Y aa' 
7 




- < 
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A similar expression is obtained for the Sine component: 

qB--qB- = ab- 

In like manner, relative physical displacements at all the nonlinear 
connecting elements are obtained, resulting in compatibility relations 
or iterating equations whose number is twice the number of nonlinear 
elements. By this process we have reduced the number of iterating 
equations from 2N (twice the number of generalized modal coordinates) to 
2M (twice the number of nonlinear connecting elements). 

Substituting the relative physical displacements (obtained by 
iteration) into the explicit expressions for the generalized 
coordinates, one then obtains the complete solution. The physical 
displacement at any point on any physical subsystem is calculated by 
simple superposition. 

The previous discussion tacitly assumed that the nonlinear connecting 
elements were planar. However, nonlinear rotor bearings deflect in a 
radial direction. This has components in both vertical and horizontal 
planes. With this regard, the displacement compatibility relations must 
be found in 2 planes, which increases the number of iterating equations 
to 4M (four times the number of nonlinear elements) - still a 
considerable reduction from two times the number of generalized 
coordinates. Note that the generalized modal coordinates include 
vertical as well as horizontal degrees of freedom. 

It will be noted that the global matrix equation of motion as well as 
the compatibility relations at nonlinear connections are derived in the 
standard way: by Lagrangian or Newtonian formulation and modal 
superposition. No recourse to Lagrangian multipliers was made to obtain 
the compatibility equations via constraint relations, as was done in 
(11). This theoretical rationalization helps establish the fundamental 
basis for the methodology. However, the results obtained in the present 
work by the standard formulation are identical to what would be derived 
with the concept of Lagrangian multipliers (11). 

3.2 Application to Turbine Engine Steady State Response 

The preceding theoretical approach was applied to the turbine engine 
problem employing the earlier formulation of TETRA. Thus the reconstruction of 
the entire engine’s dynamic response from the modal characteristics of its 
components' normal modes is made by modal synthesis. The modeling of the 
engine structure is unchanged. 

However, rather than selecting the time history of the engine response (for 
the transient case), one picks the range of frequencies where the solution is 
required in the steady state analysis. 
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Because the general rub or bearing element (nonlinear spring) produces a 
force that is a cubic in the relative displacement, the radial nature of the 
displacements means that the connecting forces couple motions in the horizontal 
and vertical planes. Thus the components of the forces and physical 
displacements at the nonlinear connections are four in number: 

(1) Vertical 

(2) Horizontal 

(3) Cosine Component 

(4) Sine Component 

This means that the number of iterating or compatibility relations is four 
times the number of nonlinear connecting elements rather than the two times in 
a system of a single plane. 

3.3 Solution of the Compatibility Relations 

The reduced algebraic equations governing displacement compatibility at 
nonlinear connecting elements is equal in number to 4-times the number of 
nonlinear connections. Principal unknowns are the relative physical 
displacements at these junctures. This system of nonlinear simultaneous 
equations are solved by iteration. Because of the doubly nonlinear character 
of the general rub element, the convergence of the iteration routine may not 
always be certain. The double nonlinearity arises from this, that the rub (or 
nonlinear) element has both a cubic term and a deadband. The deadband itself 
is a nonlinear property even when the associated spring rates are constant. 
However, the convergence problem is not really serious when deflections are 
small (in the order of bearing clearances). Only when deflections are very 
large will these problems arise. 

In the computational algorithm developed for TETRA 2, the iteration 
subroutine is written as a module which may be replaced or added to - at the 
option of the user. The program has 2 IMSL iteration modules. These are based 
on Newton's Method and the Secant Method. 

The relative physical displacements, determined by iteration, allow 
calculating all the forces in the system of equations in the generalized 
coordinates, and hence, all the modal amplitudes. The latter is performed 
simply by inversion technique. Subsequently, the physical displacements at all 
points in the complete structure are fcund by superposition. These, as well as 
the bearing loads, are found as in the original TETRA. 

The following sections are the amplification and application of the basic 
theory to the detailed analysis of the TETRA engine model. These contain the 
explicit and very detailed working relations describing the global matrix 
differential equation formation, the transformation to algebraic equations by 
harmonic balance, formulation of the iterating equations and the special 
treatment of the deadband. The latter describes the numerical procedure for 
performing the harmonic balance when rubbing is intermittent rather than 
continuous. 

In the next volume, which is also the user’s manual, sample cases are 
described along with input and output descriptions. Because TETRA 2 replaces 
the original TETRA program, with the dual capability of transient and steady 
state analyses, the full and complete set of input sheets (with descriptions 
and instructions) and output options are given. 
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4.0 Part II: Detailed Application of the Theory to Turbine Engines 


The turbine engine modeling developed in the original TETRA has been 
presented, e.g., the ordered concatenation of the structural subsystems in 
terms of their normal modes and as governed by the method of modal synthesis. 
However, for the steady state cases the equations to be solved are the 
displacement compatibility relations at the nonlinear connections, which number 
four at each of these joints. 

As discussed in Part I, the steady State response analysis requires the 
following: 

1) Formation of a global matrix differential equation in the 
generalized coordinates 

2) Transformation of these equations by the principle of harmonic 
balance with simultaneous nonlinear algebraic equations 

3) Formation of displacement compatibility relations at nonlinear 
connecting elements 

4) Iterative solution of the compatibility relations - yielding the 
relative physical displacements at the nonlinear connections 

5) Substitution of these physical displacements in the transformed 
global equations of the generalized coordinates, and calculating 
the latter by inversion 

6) Calculating the physical response of the entire engine by 
superposition 

These calculations were made in the frequency domain so that the results 
describe the forced steady state response of a turbine engine at the principal 
harmonic of the excitation frequency. 

The implementation of the methodology in the original TETRA program has 
resulted in TETRA 2, which now has the dual capability of the transient version 
and steady state. In addition, the scope of the transient version has been 
increased to contain the following enhancements initially developed for the 
steady state analysis: improved rub element (with the cubic nonlinear factor), 

structural damping capability (applies for physical connecting element types 1, 
2, 4, and 5), and new printout options. 

Because the transient analysis has been virtually untouched, this portion of 
the report is concerned only with theoretical details of the steady state 
analysis. Only in the user’s manual, (Volume 2) are the transient and steady 
state options merged together via input sheets and output. 
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The following sections document the detailed working equations and their 
implementation into TETRA 2. One should use the fundamental global matrix 
differential equations presented in Part I as a general reference in reading 
what follows, because this provides a bird’s eye view of the equations to be 
solved and the inter-relationships of the various elements. For instance, it 
should be noted that forces are initially derived in physical space and 
subsequently developed in generalized coordinate space. This is followed for 
flexible bladed disk, gyro linear connecting elements, and nonlinear 
connections. 
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4.1 Applied Forces For a Steady State Analysis 


4.1.1 Physical Unbalance Forces 

The unbalance forces are illustrated in the figure on the type K-1 input 
sheet. Only the steady state unbalance forces are discussed here, since the 
transient analysis unbalance loads were covered in reference 1. For a steady 
state run, the unbalance forces are: 


_ 2 

Fyjy = f/o) cos((ot + ({)) 


where: 

^ujz “ Unbalance load at point j in the z (vertical) direction. 

Fujy " Unbalance load at point j in the y (horizontal) direction. 

o> = rotor speed 
t = time 

U = unbalance magnitude for the unbalance load 
4> = Phase angle for the unbalance load 


Note: For each unbalance load the global point number j, unbalance 

magnitude U , and phase angle $ is input (see type K-2 namelist input 
sheet). More than one unbalance load may be inputted for a given point. 


Using trigonometric identities, we can rewrite the unbalance forces in terms of 
cos and si-n components as follows: 


Letting 


2 2 

F = C7(o sin <b cos 0 )^ + C/o) cos^sinoyt 

UJZ ^ ^ 

2 2 

F = Uoi cos d) cos (ot — Udi sin <b sin <ot 
ujy ^ ^ 

F =F^ cos (*)t + F* sin. cot 

UJZ UJZ UJZ 

F =F^ costiit + F^ sindJt 

ujy ujy ujy 


we see that the magnitudes of the cos and sin components are: 


. = Uu> sin ^ 

UJZ 


F* = Ud>^ cos4> 

ujz 


ty 

F^ = Uci cos ^ 
ujy 


F* = —U sin 4> 
ujy 
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4.1.2 Physical Pcos Forces 

For a steady state analysis run, TETRA 2 also makes provision for inputting 
applied loads of the form: 

F = P cosUat + <J)) 


where: 

^pjk = Applied forces at global point j in global direction k. 

P = Force magnitude for the load 
« = Steady state forcing frequency 
t = Time 

4> = Phase angle for the load 

Note: For each Pcos (co^ + cj)) load, the global point number j, force 

magnitude P, phase angle <J> , and global direction number k is inputted (see 

type L-2 input sheet). 

Using a trigonometric identity we have: 

F = Pcos d) cos a)^ — P sin d> sin oit 

pjk 


Letting 


F = F^ .cosoit + F* , sin(ot 

pjk pjk pjk 

we see that the magnitudes of the cos and sin components of the force are: 


F^ , = Pcosd) 

pjk 

F* , = — P sm<t> 

pjk 

4.1.3 Total Physical Applied Forces 


The total applied force for a given point and direction can be written in terms 
of the cos and sin components as follows: 

F , = F^ .cosu)t + F^ ..sinuit 
ajk ajk ajk 


where: 


^qjk = Total applied force for global point j and global direction k 

Kjk are obtained by summing the unbalance loads and the Pcos((ot + <p) 

loads for global point j and global direction k as follows: 




= Tf*-. + 

^ ujk R 


pjk 


= y F® , +F® . 

ajk ujk pjk 
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where 


^ujk * Magnitude of the cos component of an input unbalance load for global 
point j and global direction k (see section 4.1.1) 
ujk = Magnitude of the sin component of an input unbalance load for global 
_ point j and global direction k (see section 4.1.1) 

* Magnitude of the cos component of an input p cos (co< + <j>) load for 
global point j and global direction k (see section 4.1.2) 

» Magnitude of the sin component of an input Pcos{(jit + ^) load for 
global point j and global direction k (see section 4.1.2) 

4.1.4 Generalized Applied Forces 

The generalized force for generalized coordinate i may be written in terms of 
the cos and sin components as follows: 


F = cos (Jit + P* sin (of 
( ( t 


where: 

F. . Generalized force for generalized coordinate (global mode) i 

n = i i 

7 = 1*=1 
7 = U=1 

where 

ojk = Magnitude of the cos component of the total applied force for global 

g point j and global direction k (see section 4.1.3) 

= Magnitude of the sin component of the total applied force for global 

. point j and global direction k (see section 4.1.3) 

ijk Mode shape for global mode i, global point j, and global direction k 
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4.2 Linear Physical Connecting Elements 


There are six types of physical connecting elements. All of the six types can 
be used for either transient analysis or steady state analysis runs with the 
exception of the type 6 element (squeeze film damper element), which can be 
used for transient analysis runs only. The linear physical connecting elements 
include type 1 (general spring - damper element), type 2, (link - damper 
element), type 4 (engine support - links element), and type 5 (uncoupled point 
spring - damper element). This section discusses the equations used for the 
linear (types 1,2,4, and 5) physical connecting elements for a steady state 
analysis run. The equations used for the nonlinear type 3 element (rub 
element) for a steady state analysis run will be discussed in section 5. 

4.2.1 Transformation Matrix 

For steady state analysis runs, a transformation matrix [ 4> ] must be 
calculated for each of the linear (type 1,2,4 or 5) physical connecting 
elements to aid in calculating the contributions of the element to the global 
matrices. 

A sample transformation matrix is shown in figure 4-1. This figure shows the 
transformation matrix for element 3 of the demonstrator model. For this 
element, joint I (global point number 5) lies on the rotor (consisting of modal 
subsystems 1 and 2) and joint J (global point number 2) lies on the case 
(consisting of modal subsystems 7 and 8). Modal subsystem 1 (the vertical 
plane subsystem) has 5 modes (global modes 1 through 5), modal subsystem 2 has 
5 modes (global modes 6 through 10), modal subsystem 7 has 3 modes (global 
modes 11 through 13), and modal subsystem 8 has 3 modes (global modes 14 
through 16). The transformation matrix simply consists of the mode shapes 
filled in at the appropriate positions, and the remainder (the bulk of the 
matrix ) filled with zeroes. 
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ELEMENT DIRECTION 


Figure 4-1. Sample Physical Connecting Element Transformation Matrix $ 
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4.2.2 Stiffness Contributions 


The contributions of a linear (type 1, 2, 4 or 5) physical connecting element 
to the global stiffness matrix is: 


[SCI = [4>f 


where 

[SC] = Stiffness contributions matrix for the element 
[<t>r " Transpose of the element transformation matrix 
[jQ = Element stiffness matrix 

[^] = Element transformation matrix (see section 4.2.1) 

[sc] is usually a smaller matrix than the global stiffness matrix, and the 
terms of the stiffness contributions matrix [sc] must be added into the global 
stiffness matrix at the appropriate positions. 

For a more detailed description of the global stiffness matrix, see section 

4.4.3 

4.2.3 Damping Contributions 

4.2.3. 1 Non-Structural Damping 

For non-structural damping, the element damping matrix [c] is constant and is 
not a function of rotor speed or forcing frequency. The element damping matrix 
[c] for non-structural damping may be defined by input damping matrix 
definition, or computed from input damping coefficients, or computed using the 
element stiffness matrix, input element Q factor, and input element selected 
frequency. In any event, the non-structural damping contributions of a linear 
(type 1, 2, 4 or 5) physical connecting element to the global velocity matrix 
is as follows: 


[NSDC] = [<t)]^[C][4>l 


where: 

[NSDC] = Non-structural damping contribution matrix for the element 
[ 4 >]^ = Transpose of the element transformation matrix 
[C] = Element damping matrix 

[<j>] = Element transformation matrix (see section 4.2.1) 
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[NSDC] is usually a smaller matrix than the global velocity matrix, so the 
terms of the non-structural damping contributions matrix [NSDC] must be added 
into velocity matrix [Cq] at the appropriate positions. For further details, 
see section 4.4.4. 

4. 2. 3. 2 Structural Dampina 

For structural damping, the contributions of the element to the global velocity 
matrix are calculated using the element stiffness matrix [K], input element Q - 
factor, and either the independent rotor speed (if ISF = 1 on the type A input 
sheet) or the steady state forcing frequency (if ISF = 2 on the type A input 
sheet). The structural damping contributions are first collected into a 
structural damping contributions matrix [Cj] (see section 4.4.4). The 
contributions of each of the linear (type I, 2, 4, and 5) physical connecting 
elements to this structural damping contributions matrix is: 

[SDC]= ^[<t>f[K][^] 


where: 

[SDC] “ Structural damping contributions matrix for the element. 

QF » Input Q - factor for the element 
r.,r - Transpose of the element transformation matrix 
^ = Element stiffness matrix 

[$] » Element transformation matrix (see section 4.2.1) 

[SDC] is usually a smaller matrix than the structural damping contributions 
matrix [Cj], so the terms of [SDC] must be added into [C^] at the 
appropriate positions. After the structural damping contributions matrix 
[C-] is computed, it must be multiplied by either the reciprocal of the 
independent rotor speed (if ISF = 1 on the type A input sheet) or the 
reciprocal of the steady state forcing frequency (if ISF = 2 on the type A 
namelist input sheet). For more details, see section 4.4.4. 

Structural damping is a new feature that has been added to TETRA 2 , and which 
was not present in the earlier TETRA program. This feature may be used for 
either steady state or transient analysis runs. Structural damping is 
applicable for the linear (type 1, 2, 4, or 5) physical connecting elements 
only. 


4.3 Gyroscopic Elements 

4.3.1 Transformation Matrix 

For steady state analysis runs, a transformation matrix [ 4> ] must be 
calculated for each gyroscopic element to aid in calculating the contributions 
of the element to the global matrices (just as is done for the linear physical 
connecting elements). 
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A sample gyroscopic element transformation matrix is shown in figure 4-2. This 
figure shows the transformation matrix for the gyroscopic element at global 
point 4 of the demonstrator model. This gyroscopic element lies on rotor 1 
(the independent rotor). Rotor 1 consists of modal subsystems 1 and 2. Modal 
subsystem 1 (the vertical plane subsystem) has 5 modes (global modes 1 through 
5), and modal subsystem 2 also has 5 modes (global modes 6 through 10). The 
transformation matrix simply consists of the applicable mode shapes filled in 
at the appropriate positions, and the rest of the matrix filled with zeroes. 

4.3.2 Velocity Contributions 

The gyroscopic element contributions to the global velocity matrix are 
collected in a gyroscopic contributions matrix for the independent rotor [Gj] 
and a gyroscopic contributions matrix for the dependent rotor [Gq]. The 
contribution of a gyroscopic element to either the gyroscopic contributions 
matrix for the independent rotor [Gj] (if the gyroscopic element lies on the 
independent rotor) or the gyroscopic contributions matrix for the dependent 
rotor [Gq] (if the gyroscopic element lies on the dependent rotor) is: 

[GC] = [4>f [G][<i)l 


where: 

Gyroscopic contributions matrix for the gyroscopic element 
Transpose of the gyroscopic element transformation matrix 

° = Gyroscopic element matrix at point (joint) I on 

■ p ° J the rotor 

Tp = Polar mass moment of inertia at point (joint) I on the rotor 
[4)1= Gyroscopic element transformation matrix (see section 4.3.1) 

The gyroscopic contributions matrix for the independent rotor must then be 
multiplied by the independent rotor speed, and the gyroscopic contributions 
matrix for the dependent rotor must be multiplied by the dependent rotor 
speed. See section 4.4.4 for more details. 


[GC] = 

[<t>f = 

[G] = 
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GLOBAL MODE NUMBER 


o 
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4.4 Formation of the Global Matrices 


4.4.1 Matrix Equation 
The equation of motion is: 


[M]q + [K]q + [C]q = F + H 


where: 

q = generalized displacement 
[M] = global mass matrix 
[K] - global stiffness matrix 
[C] » global velocity matrix 

[F] = generalized applied (external) force vector 
[H] = generalized nonlinear (rub) element force vector 

The following sections detail how we find the [M], [K], and [C] global 
matrices. In addition, we will also show how we combine the [M], [K], and [C] 
matrices into one global solution matrix [GM]. 

4.4.2 Global Mass Matrix (Ml 


M 


M 

, . *f" 

M 




modal masses 



flexible bladed disks 


The global mass matrix consists of the diagonal modal mass terms and the 
non-diagonal terms due to the flexible bladed disks (if any). These terms are 
saved in arrays in TETRA. The global mass matrix itself is not saved in TETRA, 
but rather the terms of the global mass matrix are written directly into the 
global solution matrix in subroutine GL0B2. 

A sample global mass matrix is shown in figure 4-3. This figure shows the 
global mass matrix for a model consisting of one rotor on which two flexible 
bladed disks are located. The rotor is composed of two modal subsystems (one 
for the vertical plane and one for the horizontal plane), and each flexible 
bladed disk is also itself a modal subsystem, making a total of four modal 
subsystems. 

Note that the steady state analysis global mass has the same format as the 
flexible bladed disk mass matrix used for a transient analysis run (which is 
shown on page 70 reference 2). However, the transient analysis flexible bladed 
disk mass matrix only includes the modes for the rotor on which the flexible 
bladed disks are located and the modes for the flexible bladed disks, and is 
only found if at least one flexible bladed disk is present. The steady state 
global mass matrix, on the other hand, includes all the modes for the model, 
and is always needed, even if there are no flexible bladed disks. 
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Vertical Horitontal FBD No. 1 FBD No. 
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Figure 4-3. Sample Global Mass Matrix [M] 




The global stiffness matrix consists of the modal stiffness terms (along the 
diagonal) plus the contributions of the linear (type 1, 2, 4, or 5) physical 
connecting elements. The contributions of each linear physical connecting 
element (see section 4.2.2) are added into the correct location in the [K] 
matrix. The [K] matrix is stored separately in TETRA. Note, however, that the 
modal stiffnesses of the flexible bladed disks (if any) vary with the speed of 
the rotor on which the flexible bladed disks are located (these are the only 
stiffness terms that vary with rotor speed or forcing frequency), and hence the 
FBD modal stiffness terms are not incorporated into the [K] matrix of TETRA but 
rather written directly into the global solution matrix in subroutine GL0B2. 


4.4.4 Global Velocity Matrix 


[Cl - [Cq] + |[C^1 + Qj [G^l + 
where 

n 

( 

= the contributions of the modal damping plus the non-structural 
damping contributions of the linear (type 1, 2, 4, or 5) physical 
connecting elements. The contributions of each linear physical 
connecting element (see section 4.2.3. 1) are added into the correct 
locations in the [Cg] matrix. The [Cg] matrix is the same size 
as the global velocity matrix [C] (square matrix whose order = total 
number of modes) 

= Frequency used for the structural damping. This corresponds to the 
independent rotor speed (if ISF » 1 on type A input sheet) or to the 
steady state forcing frequency (if ISF = 2 on the type A input 
sheet) . 


[Cl= i 

= Matrix for the structural damping contributions (not including the 
multiplier) of the linear (type 1, 2, 4, and 5) physical connecting 
elements. The structural damping contributions of each linear 
physical connecting element (see section 4. 2. 3. 2) are added into the 
correct positions in the [Cj] matrix. The [C^] matrix 
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is the same size as the [C] matrix except that the rows and columns for the 
flexible bladed disks are omitted (C^ is a square matrix whose order equals 
the total number of modes in subsystems 1 through 11). 

= Independent rotor speed 


i = l 

* matrix for the contributions (not including the multiplier) 

of the gyroscopic load points on the independent rotor. The 
contributions for each gyroscopic load point i (see section 4.3.2) 
are added into the correct positions in the [Gj] matrix. The 
[Gt] matrix is a square whose order = total number of modes for the 
independent rotor. 

^£> = Dependent rotor speed 


[G^]= ^ [4>^F[Gp[<I>^] 

i = l 

= matrix for the contributions (not including the multiplier) 

of the gyroscopic load points on the dependent rotor. The 
contributions for each gyroscopic load point i (see section 4.3.2) 
are added into the correct position in the [Gq] matrix. The [Gq] 
matrix is a square matrix whose order = total number of modes for the 
dependent rotor. 

^FBD= speed of the rotor on which the flexible bladed disks are located. 

[G^^]= matrix for the contributions (not including the ^fbd 

multiplier) of the gyroscopic loading due to the flexible bladed 
disks. This does not include the terms for the flexible bladed disk 
center of gravity points, which are included in either the [Gj] or 
the [Gq] matrix. See figure 4-4 for the contents of the [GpQQ] 
matrix. 

The global velocity matrix [C] is not stored separately in TETRA, but the 
component matrices [Cq], [Cg], [Gj], and [Gq] are stored separately. 

The component matrix [SpBol is not stored, out the non -zero terms of this 
matrix are stored in arrays and scalar variables. The terms of the component 
matrices are multiplied by the appropriate multipliers and incorporated 
directly into the correct positions in the global solution matrix in subroutine 
GL0B2. 
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4.4.5 Global Solution Matrix FGM1 


Using tensor notation, our equation of motion is: 
^ij*^j ^ij^j “ *^i ^i 


where: 

= global mass matrix 
K^j = global stiffness matrix 
C^j = global velocity matrix 
= generalized displacements 
FI- = generalized applied (external) force vector 
= generalized nonlinear rub element force vector 

We can substitute: 

a . = a ■ cos (lit + b . sin (Ot 
J J 

F = F^.cos(ot + F*sinuit 
i I * 

H = cosuit + H sin(jit 

I t t 

(the expression is possible because we use harmonic averaging to get the 
nonlinear (rub) element forces, as will be shown later). 

This gives: 

— co^ Af . ( a . cos 0 )< + 6 . sin a)^ ) + K..{ a cos at + b . sin a)^ ) + aC. 1 —a . sin at + b . cos at ) 

V \ J / ij\ j j ) ‘A J j 


— cos at + sin at + H^. cos at + H* sin at 

lilt 


Now separating out the sin and cos parts this yields two equations: 

— a^M. a .cosat + K. .a cosat + aC. b cosat = F^cosat + H^cosat 

ij J tJ J iJ J I I 

— a^M .b sinat + K. .b sinat — aC a sinat = F^sinat + H^sinat 

ij J ij J iJ J I I 


Dividing through by cosat and sinat and rewriting: 

( -a^M + if .. )a . + (oC. b = 

V iJ ijJ J ij J I I 

{-a^M . + K.]b- (oC . a . = F* + H^. 

\ iJ ijJ I ij J I ‘ 
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We can write this in matrix form as follows (also noting that * 0 for i^j) 


— 


— coC 


11 


— coC 


fll 


K. 

\n 


-w M + K 

nn nn 


— coC 


1/1 


— o)C 


coC 


11 


coC 


nl 


-CO 


mC 


In 


oJC 


.K 


In 


K 


nl " 


nn nn 



«1 




^1 


a 




n 


n 

= 

n 

+ 





^1 



1 

b 

n 


F* 

n 


i 

• n 

J 



_ — 


_ 


where n = number of generalized coordinates (modes) 

We define the large 2n x 2n matrix on the left hand side to be the global 
solution matrix [GM]: 


GM\ 


l-coX + ^ll 


K 


nl 

— coC. 


11 


— (oC 


/il 


K 


1/1 


-CO M +K 

nn nn 


— coC 


In 


coC 


11 


coC 


nl 


-“^11 + ^ 11 - 


■coC < K 

nn| nl 




In 


oh: 


nn 


In 


-0)^Af +K 

nn nn 
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We see that knowing the elements of the global mass matrix 
of the global stiffness matrix the elements of the global 
matrix and the steady state forcing frequency a> , we can 
global solution matrix. 

Then we have: 


1 — -! r 

— 

• 

- 





“i 

• 

« 

• 

• 




• 

4 

• 

• 

GM 


• 

• 

9 

• 

a 

n 

• 

• 

• 

= 

n 


n 

9 

• 

• 

9 



t 

• 

« 

• 

b 


F* 


• 

• 

H* 



n 

. 


n 


■ n 


the elements 
velocity 
calculate the 
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4.5.1 General Method 


We demonstrated in the previous section that the equation of motion can be 
written in matrix form as; 



Solving for the cosine and sine components of the generalized displacements 
(a^ and respectively): 
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Expanding the T vector in the above we get: 


r ^ 


* " 






• 


• 


• 


• 

• 


a 

n 


R'= 

n 

1 

+ 



• 


• 

• 


• 


• 


• 


b 


R* 


n 


n 


» m 





^1,1^ ■■■ ^l,a + l^l ■ ■^I2n^n 

■*' ■■■ ■*'^n + U^n ■*■ ^n+u+l^l '^n+1.2n^n 

m 

^2n.l^l ••■ ■'■^2n,n^n ^2n,n + l^l ^ 2n^n^n 


( 1 ) 


where, the g^ a are the elements of the inverse of the global solution matrix 
[GM]-1. 

The g^j values (elements of the inverse of the global solution matrix) are 
found in a straight forward manner by first finding the global solution matrix^ 
[GM] as detailed previously, then inverting the global solution matrix. The 
and R* values (the R vector) can also be calculated in a straight forward 
manner' by multiplying the inverse of the global solution matrix by the applied 
(external) force vector F (also easily found), as per the definition of the R 
vector. 

If there are no nonlinear (rub) elements, or. if there is no rub because the 
dead band has not been exceeded for any rub element, then the sine a^nd cosine 
components of the generalized nonlinear (rub) element forces (the and 
respectively) are zero. For this case, the cosine and sine components of the 
generalized displacements (the a^ and b.j respectively) can be found easily 
since they are equal to the and R* values, and we’re done. However, if 
nonlinear (rub) elements are present, and if there is a rub for at least one of 
the rub elements, then we must find the h‘ and H* values, which becomes 
much more involved. ‘ ‘ 

To find the and values (if needed) iteration must be performed. This 
is done by deriving a set of iterating equations in which the rub element 
relative displacement components (the Ax’s ) are the unknowns. For each 
iteration, we calculate the rub element physical force components at joint I 
(the P’s), which are functions of the Ax’s , and which get plugged into the 
iterating equations. Then, after the iterating equations have (hopefully) 
converged to a solution, we take the final values for the rub element physical 
force components at joint I (the P’s), and calculate the nonlinear (rub) 
element generalized force components (the and ff* values) (the H vector). 
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The rest is simple. The T vector is then found by multiplying the inverse of 
the global solution matrix by the H vector as per the definition of the T 
vector. The R vector and the T vectors are then added together to find the 
cosine and sine components of the generalized displacements (the a^ and b^ 
respectively) . 

The following sections detail the derivation of the iterating equations and the 
derivation of the rub element physical forces at joint I. 

4.5.2 Derivation of the Iterating Equations 

We can write: 


AJf, = AX‘, cos (jit + AJTt sin (lit 
u k^v k^v 


^ cos (iit + AX* ^ sin (jit 


where ^^k,u - relative displacement for the k’th nonlinear (rub) element 

in the vertical direction 

^^k,H * relative displacement for the k’th nonlinear (rub) element 
in the horizontal direction 

For each nonlinear (rub) element k, we have four equation as follows: 

( 2 ) 


( 3 ) 


( 4 ) 


( 5 ) 


AXt =Xt ,-Xl ,= ya((j). ,-4>t , 

j = l 

^k,n ~ K,h,i ~ K,h,j ~ — “ '^*, 1 ,//, 

( = I 

fl y 

AX* =X* ,-X* = y 6 (4>^ ,) 

k,U kyV^l kyV,j ^kjjl^VyJ ] 

i-\ 

- / 

4>. w , - <J> 

*,// k,H,l k,H,J ^ 


where 


^ or s 
k, vor H, lor J 

or H,IorJ 


bi 


cosine(c) or sin (s) component of the vertical (V) or 
horizontal (H) displacement for joint I (I) or joint J 
(J) for the k’th nonlinear rub element 
Mode shape for k’th nonlinear (rub) element, mode i, 
vertical (V) or horizontal (H), joint I (I) or joint 
J (J) 

cos component of the generalized displacement for 
mode i 

sin component of the generalized displacement for 
mode i 
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We now define: 


^k,i,o,I ^k,i,v,J 

Substituting this into equations 2 through 5 we have: 

n 

£iXl = y a. A4). . 

t = l 
n 

^k,H= 


i = l 

n 


Aj:f = y 6. A4). , 

kyV ^ I 
|S»l 

n 

^;,«= 

i = l 

Further, we can define 

n 


Fl = 

k,v 

y a. A<^. . - , = 0 

i=l 

n 

(6) 

= 

k,H 

n 

(7) 

C 

II 

y 6 Ad). . -AXl =0 

^ j ^k,i,o k,o 

i = l 
n 

(8) 

II 

i = l 

(9) 


where = the iterating equation for the k’th rub element, cosine (c) 

or sin (s) component in the vertical (V) or horizontal (H) 
direction 

The unknowns in the iterating equations are the relative displacement components 
iK,*-’ K,h) • During iteration the values of the ‘'Xu- ‘X,h ) 

express! ohs approach 0 as the unknowns approach their true value. The number 
of iterating equations equals 4 times the number of nonlinear (rub) elements. 
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From equation 1 we see that: 

j=i j=i 


J=l 


j=l 


Plugging the a^ and b^- expressions into the iterating equations (equations 
6 through 9) we get: 


II 

M 

R‘ + y g . Ji‘ + g ^ n‘^ 

• '"'I'm..- 

= 

k,u 

1=1 

7=1 7=1 



" / 

n n 



= 2 

R'+Yg ir+\'g 

1 1.7 7 °i,n+j Jj 


- AX^ 

k,H 

(=1 

7=1 7=1 



n 

, n n 

\ 


= 2( 


)a«i>, . 

f k,l^u 


j=i j=i 

n n 




.=1 j=i 


7=i 


Now defining: 


> R‘ 

I kyV 

( = l 

~ ^k,H 

i = l 

fl 

y Act), = S* 

( ^k,i,u k,v 

1 = 1 


)=i 


We can write: 


( 10 ) 

(11) 


''♦..i,, I 1 I AX', „ = o 
+ 1 i «,,//; + i A4. , y e H‘ - AX' = 0 

i-l ^=1 J_^ > J J 

n. = i + 1 A4., , _ y ,, //• - AX^, = 0 ,12) 

^'m/ = + i. ^'^k^ui i + + 1 i 5,,^, „ //; - AX^ = 0 (13) 

1 = 1 J=l , = I _ 1 ’ j J x.ll 
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To get the generalized nonlinear force (the H’s) we add the contributions of 
joint I (subscript I) and joint J (subscript J) in both the vertical (subscript 
V) and horizontal (subscript H) directions and sum over the nonlinear elements: 


//*'=> , + //%,, + Fi‘ , + , 

J — \ J,0,I J,n,l J,V,J j,H,J 

e=i 


, + //* „, + //* , + //" „ , 
J ■— \ J,U,I J,n,l J,U,J J,H,J 
t=l 


To get each H contribution we multiply the mode shape *1* by the nonlinear 
connecting element force P. 




e=\ 



Noting that: 



pC _ _ uc 

^e.u.J i,u,i 

P‘ = 

e,H,i 

p* = _ p* 

^e,v.j (,o,i ’ 

ps = 

' e,n,J 

-P'* 

e,H,i 


we get: 


e=i 


Defining: 


A<I>, = r - , 

“jfV c,7»u,/ c,7, 


^e.j,H,j 
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we get: 


(14) 


^ / 

, A(I>, + AO 1 

W* = V ( P* AO + P* AO 

j ^ ^ ^(,Hj ej,H 


(15) 


Substituting these expressions for //‘'and //'’ into equations 10 through 13 we 
get: ^ ' 


I) ~ ^ a4>. 8 'y 

*1 u •*— ~k,i,u ° i,j *“ 

^=1 r=i 


f-i 

n 


l.y./ ^ 


+ AO, . 'y 8 ^ ( P“ Ad) + P'* A(b 

— ^i,n+j — \ e,u,i ^e,j,u e,nj , h 

1 = 1 j = i f=i 


-^X^. = 0 

», II 


"In- ^llJ + 1 + ''^.//.Z „ 

J = 1 ^=1 '-v.” 


1 = 1 


1 = 1 


+ I I I AO + PI AO 

^ = 1 f=l 


-‘^U, = « 


S ‘•'^l.j.H 






-AZ? =0 

« , 9 


^’a,h + i. ^4>t X ( ^e,u,/ ^4>^,v.y ■*■ » ) 

1 = 1 ; = 1 e=l ' 


n n L 

N A(b V V 


1 = 1 


J=l 


f=l 


+ 2. + ^u/,/ 




-^1// = 0 
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Rewriting 


n 

n 

L 

n 

n L 



u ~ “1,7 

y a 4 ). 

•^' c, 7 ,i) 

K , + y ‘^ 4 >i. 

U,/ — ^*,1 

,11 —^1,7^ ^‘^e. 7 .// ^t,n,i 

( 16 ) 

i = l 

7=1 

f=l 

i = l 

7=1 e=i 


n 

n 

L 

n 

a L 



+ N' Ad), >5 y Ad), . p", , + A<t), y g ^ y a 4). ,, p" ,, , 

^ ~ ® I,n+J -A- ^e,j,v e,u,l — ^k,i,u — ° i,n+j ^ — t,J>" c.HJ 

j=\ e=i i=i 7=1 e=i 


1 = 1 


-AXt = 0 

X, U 


j=l ?=1 i=l 7=1 ^=1 


i = l 


X ^ ^i ,«+7 ^ ^ ^i,a+ 7 ^ ^t,HJ 

J=l e = l i = l 7=1 ^=1 


i = l 




n L n n L 

= Si + y Ad), y g y a*, pi , + A(1), ^ g ^ y ^4), „pi „ , 

j=i e=\ i=i 7=1 ^=1 


1 = 1 


(18) 


n n L n n L 

+ X ^i,n+7 ~ ^‘^f,7.i- ^?.u./ ■*■ ^ ^ ^n + i,n+j ^ ^‘^e. 7 , » 

1 = 1 j=i e=i 1=1 7 = 1 ^=1 

-AXl =0 

ik, V 


n L n n L 

K,II ~ X ^‘^’*,1,// ^ ^ri + 1,7 ^ X ^\i,H — ^n + 1,7 ^ ^ t, H ,I 

J=\ e=l i=l 7 = 1 ^=1 


i = l 


n n L n n L 

^ S ^n + i,n+j ^ ^‘^?, 7 .u ''' — ^yn + i,n +7 ^ ^‘^e. 7 .« 


i=l 


’ n + i, n+j 

7=1 ?=1 


1 = 1 


' n + i, n+7 

7=1 f=l 


(19) 


-^u/ = o 
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We define the following parameters: 


A4>, . yg A<i>. . 

i=l j=l 

n n 

i=l j=l 

n n 

A^ = y A<b B ■ Ad). 

k,o,e — ^k,i,u ° n + i,j 

i=l J=1 

n n 

(=1 J=1 

fz n 

= y A<p, . y g A<p. „ 

kyUyt — ^kylyU ^lyj ^Zyjyli 

i = l 7=1 

n n 

i=l 7=1 

n n 

K,o,e = I^‘1>*.,,.S5„+..7 

i=l J=l 

n n 

^k,HJ ~ ^ X ^n + 1,7 

.=1 7=1 

fi n 

(f, , = y A({> . yg , A<i>. 

/fe, y, C — ^ kylyO ^ ^ lyTl^-J tyJyV 

t=l J = \ 

n 


1 = 1 
n 


7 = 1 
n 


1=1 7=1 

n n 

y 


^ u e ~ A<j>, yg Ad). 

— '. k,i,u ^ ^ n + i,n+j ^e,j,u 

c* 


'k,H,e 


( = 1 
n 


; = I 

n 


^ rt + 1 , n+j 


^*, V, e X ^^k, i,u^ ^i,n+j H 

i = I 7=1 

n n 

^k,H,e~ ^ ^^k,i,H 'E. ^i,n+j ^‘*’^,7.// 
,=1 7=1 

n n 

u,e^ E ^^k. i.vEs„+i^„+j H 

.=1 7=1 

n n 

^k,H,t ~ E ^‘^k,i.U E ^n-^i,n+j 

1=1 7=1 
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( 20 ) 


Substituting these parameters into equations 16 through 19 we get: 

L 




+ cf. . Pi . , + Di 


-AX, =0 


’ ’ f=i 

?=i 

— o* 4. ( 4 pc + H'* + D“ P“ ) — AX" = 0 

u~ ^k,v^ — \^k,u,e ^e,u,i ^ °k,a,t ^ ^kuj ^ e,u,i k,u,e e.H,i ) k,u 

f=i 


n 


( 21 ) 

( 22 ) 

(23) 


t=i 


These are the iterating equations. The kji • ^k,v’ ^kji 

functions are calculated in subroutine FCN. The actual iteration is performed 
either by IMSL subroutine ZSCNT or IMSL subroutine ZSPOW, depending on which of 
these routines the user specified via the nonlinear routine option IROUT (see 
input sheet 1-2) . 

The S, A, B, C, and D parameters in the iterating equations are known 
quantities and are calculated based on the proceeding definitions of these 
quantities. Because the S, A, B, C, and D parameters do not change from 
iteration to iteration, for efficiency purposes they are calculated only once 
for each forced frequency, prior to iterating for the forced frequency. 

The are the unknowns in the iterating equations. The P’s in the 

iterating equations (the rub element forces at joint I) are functions of 
the AX’s . See section 5 for a discussion of the method and the 

equations used to find the P’s. 

4.5.3 Equations for the Maximum and Minimum Relative Displacement 
Magnitudes 

The maximum and minimum relative displacement magnitudes for each rub element 
are needed to determine (through comparison with the dead band) if there is a 
rub for the given rub element and, if so, whether the rub is continual or 
intermittent. 

As noted previously: 

AX, = AX^ cos to/ + AX" „ sin uil 
H = H H 
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where 

^ku ~ 'f'elative displacement for the k’th nonlinear (rub) element in the 
vertical direction 

^k,H * relative displacement for the k’th nonlinear (rub) element in the 
horizontal direction 

We next define: 


AX"' = 

V 

AX"‘ = 

k, H 




4 ) = tan 


AX* 

— I k,u 

AX^ 

k, V 


- 


= tan 




AX^; 


k,H 





max 

min 


= Maximum relative displacement magnitude for the 
element 

= Minimum relative displacement magnitude for the 
element 


k’th nonlinear (rub) 
k’th nonlinear (rub) 


It can be shown (through the derivation is long-winded, so we skip it here): 
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4.5.4 Determining If a Rub Is Present (That is. if Iteration is Needed) 


Iteration is needed to find the AX’s (rub element relative displacement 
components) only if there is a rub for a least one of the rub elements. The 
procedure used to determine if a rub is present is as follows. First, the S 
parameters are found (subroutine PARAMl). Then, in subroutine CHECK, we find 
the AX’s that would result assuming there were no rub element forces 
(that is, if the P forces were all equal to 0). From equations 20 through 23, 


these are: 


( 26 ) 

1 

1 



( 27 ) 1 




( 28 ) 1 

f (assuming no 
I rub element forces) 



( 29 ) J 



Then, using these ax’s (calculated assuming no rub element forces), we 
calculate the maximum and minimum relative displacement magnitude for each rub 
element that would result assuming no rub element forces from equation 24 and 
25. By comparing the maximum relative displacement magnitude with the dead 
band for each rub element, we can check (as we do in subroutine CHECK) if a rub 
is present for any each rub element. 

If this check reveals that there is no rub for any of the rub elements, then we 
know that the assumption of no rub element forces was correct. If this is the 
case, then the aX’s obtained from equations 26 through 29 are in fact the 
correct aX’s , and no iteration is needed. On the other hand, if this check 
reveals that there is a rub for at least one of the rub elements, then there 
are rub element forces, and equations 26 through 29 do not yield the 
correct AX’s . For this case, the AX’s must be found via iteration 
using equations 20 through 23. 

4.5.5 Finding the Initial Guess for the Iteration 

Assuming that we have found that iteration is needed (see section 4.5.4), then 
we must find the initial guess for the iteration. This is done in subroutine 
check. The initial guess is determined from the following rules: 

1. If it is the very first solution (forced frequency), and the user inputted 
an initial guess via the GUESS input variable (see namelist input sheet 
1-3), then the initial guess inputted by the user is used. 

2. If it is the very first solution (forced frequency), and the user did not 
input an initial guess, then the initial guess for the rub element 
relative displacement components are those that would result if there were 
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no rub element forces (obtained from equations 26 through 29). 

3. If it is not the first solution, but if rotor speed is considered for the 
run, there are no unbalance forces, and if it is the beginning forcing 
frequency for the current rotor speed, then the initial guess for the rub 
element relative displacement components are those that would result if 
there were no rub element forces (obtained from equations 26 through 29). 

4. If rules 1, 2, and 3 do not apply, and if the previous frequency had a rub 
for at least one rub element, then the rub element relative displacements 
(found via iteration) for the previous forcing frequency are used as the 
initial guess for the iteration. 

5. If rules 1, 2, and 3 do not apply, and if the previous forcing frequency 
did not have a rub for any of the rub elements, then the initial guess 
for the rub element relative displacement components are those that would 
result if there were no rub element forces (obtained from equations 26 
through 29). 

4.5.6 Solving the Iterating Equations 

Assuming that we found that iteration was needed (see section 4.5.4), the first 
step performed is to find the initial guess for the rub element relative 
displacement components (see section 4.5.5). Next, we proceed to find the A, 

B, C, and D parameters (subroutine PARAM 2) that appear in the iterating 
equations. Then a whole bunch of subroutines are called (if needed) to solve 
the iterating equations (subroutines SOLVE, NLFORP, BACKS, FCN, plus IMSL 
subroutine ZSCNT or ZSPOW and several other IMSL subroutines or function 
subprograms which are called by ZSCNT or ZSPOW). If there is no rub for any 
rub element (so that iteration is not needed), all of these subroutines are 
skipped. 

The user specifies which of two IMSL subroutines (ZSCNT or ZSPOW) are used to 
solve the iterating equations via the nonlinear routine option IROUT (See 
namelist input sheet 1-2). By this means, the user chooses the iteration 
method to be used. 

At each iteration, the rub element physical forces at joint I (the P forces) 
must be found so that they can be plugged into the iterating equations. These 
forces are found in subroutine NLFORP. The next section details how we find 
these forces. 
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5.0 General Nonlinear Rub Element 


5.1 Physical Force equations at Joint I 

We now turn our attention to the equations for the rub element physical forces 
at joint I. These equations are nonlinear, both because of a cubic term and 
because of the dead band. 


We define: 


= Physical forces for the 1^^ rub element in the vertical direction 
at joint I 


= Physical force for the 1^^ 
at joint I 


rub element in the horizontal direction 



= Relative displacement magnitude for the 
fth rub element. 


AX = Relative displacement in the vertical direction for the 1^*^ rub 

element (that is, the vertical displacement at joint I minus that at 
joint J) 


* Relative displacement in the horizontal direction for the 1^^ rub 
element (that is, the horizontal displacement at joint I minus that 
at joint J) 

AX^ y = Relative velocity in the vertical direction for the 1^*^ rub element 

’ th 

Relative velocity in the horizontal direction for the rub 
element 


^0 = Dead band (this equals input variable DBAND on the type F input 
sheet) 

= Linear radial spring constant factor for the 1^^ rub element (this 
equals input variable SK on the type F input sheet) 

= Nonlinear radial spring constant factor for the 1^^ rub element 
(this equals input variable AK on the type F input sheet) 

Cf = Damping coefficient for the 1^*^ rub element (this equals input 
variable CC on the type F input sheet) 


We can write: 



for 


^ ^0 - 



for 

j 

> ^ 0 ^ 






p 

^-K.AX. I - 




e.Hj 

e e,/i\ 


t\ 


A. 


- e. 


- e. 


- C. AX^ 

c c, t> 




(30) 


( 31 ) 
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5.2 Harmonic Averaging 

In order to find the ■ ^e.u,/ ■ and ^e.n,i values 

to plug into iterating equations 20 through 23, we approximate the rub element 
physical forces at joint I with the expressions: 

^eju = ^ 


These expressions are not exact because, due to the nonlinearity and complexity 
of equations 30 and 31, higher order terms involving cos 2 w/, sm 2 io/, 
cos3o>/, sin3a)/,etc. would al SO be present. However, these higher order terms are 
neglected and only the first harmonic terms retained. Our problem then boils 
down to finding the values of - and , • This 

is accomplished using the method of narmohi’c averaging. 

Given a function f(t), the method of harmonic averaging involves integrating 
over a cycle as follows; 

f 2ii/(j 

for the cos part: fWcosuitdt 

JO 

r 2n/(i> 

for the sin part: f{t)siru^idt 

Jo 


In order to perform the harmonic averaging, the following integral 
transformations will be helpful. The derivation of these integral 
transformations is straight forward but is fairly lengthy so is not included 
here. 

For the expression: 

</(/) = a cos (A)t + b sin (oC 


We have: 


2n/o> 


q (t)cos itildl = 

0 

n 

— a 
0) 

(34) 

r 2li/cj 

q U) sin o)l dt = 

- b 

(35) 

j 0 

a> 

r 2n/cj 

• 

q (t)cos(x>tdt = 

n6 

(36) 

0 



r 2ii/a> 

q U) sin uit dt = 

— na 

(37) 

>0 

r 2ii/oj 

3iia^ a“+ 6" j 

(38) 

I q^ U) COS (tit dt = 

0 

4o) 

r 2ii/o) 

3nb{a^+b'^^ 

(39) 

q U) Sin (tit dt = 

J 0 

4w 
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For the pair of expressions: 

q^(t) = a cos Cx>t+ b sin at 
(t) = c cosat + d sin at 

we have: 


r2n/u 


q^it) c^{t)cosatdt = 


a{dF + 3c^) 4- 2bcd 


4o) 


' 2n/a> 

qM) qlit)sinatdt = 

Jo ^ ^ 


b(f? + ZdF) + 2acd 
Aa 


(40) 


(41) 


Note that the quantities, 

which are found via harmonic averaging {except when there is^no rub), ^ are ^ 

functions of the relative displacement components 

The AX’s are the unknowns in the iterating equations^, and each iteration 
provides a new guess for the AX’s . Thus, the ^e,u,r 

and values change for each iteration, so the harmonic averaging must 

be done for each iteration. 

If the equations are simple enough, closed form equations may be derived to 
perform the harmonic averaging (as is done for a continual rub with dead band 
equal to 0). If the equations are too complex to solve in closed form, 
numerical integration (using Simpson’s rule) is performed instead (as is done 
for a continual rub with dead band not equal to 0 and for an intermittent 
rub). The advantage of using closed form equations, if possible, is that the 
closed form equations are more efficient and more exact than the numerical 
integration. 

Further details of how the harmonic averaging is accomplished is contained in 
sections 5.3.2, 5.3.3, and 5.3.4. 
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5.3 Four Possible Rub Categories 


Four possible rub categories are recognized by TETRA 2. These categories are 
no rub, continual rub with dead band equal to zero, continual rub with dead 
band greater than zero, and intermittent rub. The categories are illustrated 
in figure 5-1. 

The TETRA 2 program determines which of the four categories apply by comparing 
the maximum and minimum rub element relative displacements (calculated using 
equations 24 and 25) and the rub element dead band (input variable DBAND on the 
type F input sheet). Note that this determination must be made for each rub 
element and at each iteration (since the calculated maximum and minimum 
relative displacement magnitudes change from iteration and iteration). 

Different logic is ^sed to calculate the rub element physical force 
components ^ t.H.i depending on which 

category applies. The following sections detail the equations used for each of 
the four rub categories. 

5.3.1 No Rub 


If it is determined that a given rub element has no rub (see section 5.3), then 
the rub element physical forces must be zero. 

Hence, the program sets: 


t.H.I 


= 0 
= 0 


P 


a 

e.H,i 


= 0 


and we’re done. 
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Clearance Circle 


Orbital Ellipse 


CONTINUAL RUB 

WITH 



Orbital Ellipse 


CONTINUAL RUB 

WITH 



Orbital Ellipse 


Clearance Circle 


INTERMITTENT RUB 

(shaded areas indicate 
where rub occurs) 



Orbital Ellipse 


Cl earance Ci rcl e 


Nomencl ature: 

^0 = Dead Band 

j = Relative displacement magnitude for the k'th rub element. 

= Relative displacement of the k'th rub element in the horizontal 
direction 

= Relative displacement of the k'th rub element in the vertical 
direction 


FIGURE 5-L. FOUR POSSIBLE RUB CATEGORIES 
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5.3.2 Continual Rub with Dead Band Equal to Zero 


If it is determined that we have a continual rub but the dead band equals 0 
(see section 5.3), harmonic averaging is used to find ^e,v,i > 

and • However, the equations simplify enough when the dead band 

equals o’t’hat we can derive closed form equations for P ^ „ j . P( h,j’ ’ 

and H I » avoiding the less efficient process of integrating 

numerically’. ’ To find the desired closed form equations, we proceed as follows 


Setting equations 30 and 31 we get: 


(42) 


Making use of the coscot integral transformations (equations 34, 36, 38, and 
40) equation 42 becomes: 

ay f,V,I Cl) * *•*' 4(0 ^ ^ V “ 




+2ax;,„ax5 „ ^ x •„ 


-nC.AX* 

c c, w 

Making use of the cosat integral transformations (equations 34, 36, 38 and 
40) equation 43 becomes: 




n 


4 » 




(ax; „)H 


+ 2^* „AX^ AX* 

c, n c, y c, y 


Making use of the siiuat integral transformations (equations 35, 37, 39 and 
41) equation 42 becomes: 


-P* =--X, AX* X p AX! [(ax^ y + 

e,v,i e t,u <,u[V t,oj \ t,u 


4(0 


AX* 

f, y 


(ax; + 3(ax; „)=' 




+ nC,AX;^, 
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Making use of the sintut integral transformations {equations 35, 37, 39 and 
41) equation 43 becomes: 




AX 




+ AX 








(ax] y + 3 {aX‘ ^ 


+ 2ijr, „ AX 


8 

^.0 


+ nC,AX»,^„ 


Rewriting the last four equations we get: 

-K,AX',„-<aC,AX>,. 

K.H.,= -I [^^«K'^^0' + 3('^.«)'+3(aX;„)^+(ax;.)^1 +2AX',«AX'^^ 

3(aT, „)2 + 3(ax; +3(ax; „)^ +(^e,»)i + 


^(,0,1 I 


-K^AX]^-<^,AX^,„ 




AX^ 

e,H 


3(aX; + 3(aX',,„)^ +3(aX; +(aX', ,)^| + 2AX‘, „ AX^ , AX', „ 




5.3.3 Continual Rub with Dead Band Not Equal to Zero 


If it is determined that we have a continual rub and the dead band does 
equal 0^(see section 5.3), harmonic averaging is used to find Pg^u,i ’ ^\h,i ■ 
and • The equations for this category of rub are too ’complex to 

solve in closed form, so we must integrate numerically to find 
P*. . . and Pt 


t,U,I 


not 


pC 

e.H,i 


pS 

e,v,i 


e,o,i 




Using the integral transformation of equation 34 on the expression given in 
equation 32 we get: 

^ r 2n /co 

p® , = - p, f cos (Of dt 

e,o,I n o 
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Using the integral transformation of equation 34 on the expression given in 
equation 33 we get; 




{■> 

n 


2n/u 

0 


P, coscot dt 

V j ii t i 


Using the integral transformation of equation 35 on the expression given in 
equation 32 we get: 


2n/(o 

P, , sinuit dt 

Using the integral transformation of equation 35 on the expression given in 
equation 33 we get: 



« f 

= n j 

2n/u 

P, „ , sin at dt 

0 




Substituting 'P = <ot , the 

preceding 

four equations can be rewritten 

=i r 

n Jc 

2n 

P, , cos Vci'P 
, t.uJ 

(44) 



t,H,I n 


(45) 



P' =M 

2a 

P, , sin '¥dW 

0 

(46) 



P* = - 

r2n 

P,„ ,sin 'Vd'V 

0 

(47) 



P cos 'P P cos W 

For a continual rub, we note that the products 

^<, 0,1 and repeat themselves every 180 . 

Thus, we need only integrate between 0 and n and double the results as 
follows: 

P*= 

n , 


(48) 





(49) 



P* 

(" P. , sin VciV 

p n 

(50) 



p* = - 

IS 

P, „ , sin Vci'P 
Jo 

(51) 
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To perform the numerical integration, the integrals are divided into 10 
subdivisions, with 'P varying from 0 through n in steps of n/lO.^ 

For each value of , the parameters v’ u’ h 

are calculated using the equations: , . ,u 

AX, = AX^ cos<af+AX* sin (at 
y c, u 

AX^ ^ = AX^ ^ cos (Of + AXj jj sin (at 
AX, = (oAX* (MS(at — (a^Xi sin(at 

c, y c, y y 

AX^ AX* jj(x>s(at — (a jj sin (at 


Substituting V = oj< and rewriting: 

AX, =AX^ cos'P + AX* sinV (52) 

r, y c, u c, y 

AX^ ^ = AX^ ^ cos + AX; ^ sift >P (53) 

AX^ „ = AXJ. u 

AX^ ^ = < 0 ^ AXj ^ cosV - AX^ ^ sift'p) (55) 


where « = forcing freouency, t = time^ 'P=cof , and the relative 
displacement components AX^ ^, ^t,H’ ^t,v AX*^ are known because they are 
the guesses for the current iteration. 

P P 

The and e,H values are then calculated for each value of *P using 

equations 30 and 31. 

Finally, knowing the ^t,v and e,H values for each value of >P between 0 
and n (in steps of n/io ), the values of the integrals in equations 48 
through 51 are obtained using Simpson’s rule. 

5.3.4 Intermittent Rub 

If it is determined that we have an intermittent rub^fsee section 5.3), 
harmonic averaging is used to find ^e,u,r 

Again, the equations for this category of rub ’are too’complex to’ solve in 
closed form, so we must integrate numerically to find ^t,v,r 


|S 

t,H,I ' 
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As for a continual rub with dead band equal to 0, the following equations are 
applicable (see equations 44-47 of section 5.3.3): 



P 


c 

e,H,t 


P 


s 

f, oj 


p 


s 


1 

n 

1 

n 

1 

n 


r 2n 


P, .cosVdV 

0 

• 2n 

P. „ ,cos'Prf¥ 

0 

f2n 


0 



P, .sin'PdV 

f, U, I 


For an intermittent rub, the orbital ellipse intersects the clearance circle at 
four points (points A, B, C, and D in figure 5-2a). Rub element forces are 
present only for the two portions of the orbital ellipse that rub (between 
points A and B and between points C and D in figure 5-2a). Furthermore, the 
products P f ^ j sinW , and PfjfjSinW are the same for 

the two portions that rub’ (see sample versus 'P plot in 


figure 5-2b). Hence, we need only 

between points A and B) and 

double 

become : 

2 

rVB 


P‘ 
t, v,I 

n 

P, .cosWd^ 

(56) 
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■'VB 


P^ = 

i,H,I 

n , 

P, „ .cos^d^ 

'VA 

(57) 


2 1 

■VB 


p3 _ 

e, oj 


P, .sinVdW 

VA 

(58) 


0 

rWB 


P* = 


P. ^ rSinVdV 

(59) 

t,H,I 

n . 

-PA 



Thus, the preceeding equations 


where is the angle for point A and 

point B. 


is the angle for 
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Figure 

(Shaded 


-2a. Orbital Ellipse and Clearance Circle 
areas indicate where rub occurs) 



Figure 5-2b. Sample „ / 


cos'V Versus 'V Plot 



To perform the numerical integration, the integrals are divided into 10 
subdivisions, with varying from to in steps of ( '^b “ ^ 

The procedure is the same as for the continual rub with dead band not equal to 
0. That is, first the^; „. ^e,H‘ ^f^.and ^ parameters are found for 
each value using equations 52 through 55. Then, the „ and ^ 

values are calculated for each value of using equations 30’and 31.' 

Finally, knowing the „ and P^ jj values for each value of ^ , the 

values of the integrals in equations 56 through 59 are calculated using 
Simpson’s rule. 

So far we have not covered how to find the '*'a angles over which the 

integration is performed. This is explained in section 5.4. 
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5.4 Equations For the Intersection of the Orbital Ellipse and the 
Clearance Circle 


For an intermittent rub, we must find the equations for the intersection of the 
orbital ellipse with the clearance circle. This is needed because we must 
numerically integrate over one of the two rub areas, and we must know the 
beginning angle and the ending angle ’Pg for the integration (see 

section 5.3.4). We proceed as follows: 

The equations for the orbital ellipse of a rub element are: 

^X =AX^cos'P+ AX* sin 'P (60) 

a u u 

AX„ = AXt, cos »? + AX*„ sin 'F (61) 

MM n 


where = co/ 


Defining: 


AX” = \/(aX^]^ + (aX*)^ 


4 > = tan 

u 

= ton ~ ^ 


AX* 
0 

AX* 


Equations 60 and 61 can then be written as: 


AX =AX”cos(‘y - <t> ] (62) 

u u \ u/ 

AX^ = AX” 0)S ( W - (63) 
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For the four intersection points of the orbital ellipse with the clearance 
circle (whose radius = the dead band) (see figure 5-2) we can write: 

Plugging equations 62 and 63 into equation 64 we get: 
cos* (v-«)+(AX")*c<»*(v-<t.„)=(e,)* 

Dividing through by we get: 



Using a trigometric identity this becomes: 



Rearrangi ng: 



Using another trigometric identity this becomes: 

AX'" X / X , , / AX'" 


cos 2V cos 2«I>„ + sin 2V sin 2<I>„ ' “ ^ 

rt ti i 


{-^)* +(-7^)* + +5in 2V sin 2<t> 

\ \ J \ 2\ J \ 

,i(^) 

2' «0 ' 

Rearrangi ng: 

/ AX'" X / AX?; 

l-^fcos2<P +(— ^ 

/AX'"\ /aA„x 


COS 2<I) 


1 
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H 

J 

/AX'” 

' j 
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COS 2V 


sin2W = 1 - - 
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AX'" 
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AX 
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a: 3 



This can be rewritten: 


Acos2'¥ + Bsin2^> = C 


( 65 ) 


where 
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Defining: 


-1 


B 


« = tan — 
A 


Equation 65 can be rewritten; 


Va^ + B^ :os(^2W - a ) = C 


from which we get: 


2'P _ a = cos 


-I 
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Defining: 

0 = cos 
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( 66 ) 


where: 


0 < 0 s n 



We can then list four possible values of the function: 


-i 

cos 


, g ] 


as follows: 
C 


cos 


-1 


-if 

cos I 


cos 


-1 


cos 


-1 


C 

^/FTF 

c 

Va^TW 

c 


= 0 


= 2n - 6 


= 2n + 0 


= 4n — 0 


Va^TW 

Plugging these four possible values into equation 66 we get: 
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1 2 2 

0 oc 
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2 2 2 

0 oc 

'P, = n + - + — 

3 2 2 

0 

ip = 2n - - + - 

4 2 2 


The preceding equations give the 'P angles for the four intersection points 
of the orbital ellipse with clearance circle. As discussed in section 5.3.4, 
we numerically integrate over one of the two rub areas, then double the 
result. However, given the four intersection points, there is the possibility 
that the rub areas are between points 1 and 2 and between points 3 and 4, and 
there is also the possibility that the rub areas are between points 2 and 3 and 
between points 4 and 1 (see figure 5-3). To find which of these two 
possibilities apply, we first calculate the angle at point 1.5, which is 

half way between points 1 and 2: 


'P 


1.5 


V + *P 
1 2 

2 


a 

" 2 
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V 



Figure 5-3. Orbital Ellipse Possibilities For 
An Intermittent Rub 
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We then calculate the radius at point 1.5 as follows: 



If the radius at point 1.5 is greater than the dead band, then the rub areas 
must be between points 1 and 2 and between points 3 and 4. If this is the 
case, we numerically integrate between points; 


as detailed in section 5.3.4. 

However, if the radius at point 1.5 is less than the dead band, then the rub 
areas must be between points 2 and 3 and between points 4 and 1. If this is 
the case, we numerically integrate between points: 


and = ^3 


as detailed in section 5.3.4. 
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5.5 Equations For the Maximum and Minimum Rub Element Harmonically 
Averaged Force Magnitude 


As noted in section 5.2, the rub element physical forces are expressed as: 

P. J = P\ r, 

I 


where P" p* and P* 

averaging. 

We next define: 


are found using harmonic 




pm V + f P* 


_1 #.»./ 

<t> = tan 

w pc 

t,u,I 


pS 

* , -1 

t,H,I 


$ = (j> - 4>„ 

V n 


Fmax = Maximum rub element harmonically averaged force magnitude 
Fmin = Minimum rub element harmonically averaged force magnitude 
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It can be shown (though the derivation is lengthy, so we skip it here): 


F 

max 


F . 

min 







Note that these equations are similar in form to the equations for the maximum 
and minimum relative displacement magnitudes given in section 4.5.3. 
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6.0 Generalized Displacements and Generalized Velocities 

6.1 Generalized Displacements 

The generalized displacements can be written in terms of the cos and sin 
components as follows: 

Zi^=Z^j^cos(at + Z*j^sin(,it (67) 

where: 


^ 1 , = generalized displacement for global mode k 

= amplitude of the cos component of the generalized displacement 
mode k 

2* = amplitude of the sin component of the generalized displacement 
mode k 

w = steady state forcing frequency 
^ = time 

The amplitude of the cos and sin component of the generalized displacement for 
each mode is calculated using the method outlined in section 4.5.1. The 
generalized displacements can also be expressed in terms of the magnitude and 
phase angle as follows: 



For each forcing frequency, the generalized displacements are found as outlined 
above and in section 4.5.1. The generalized displacements are then used to 
calculate the physical quantities (physical displacements, physical velocities, 
physical connecting element forces, etc.) as outlined in section 7 which 
follows. It is normally the physical quantities that the user is interested in 
rather than the generalized values. For this reason, printout of the 
generalized displacements is omitted if the user requests the short or the 
standard form of the output (see printout option lOUT on the type A input 
sheet). Printout of the generalized displacements (cos components, sin 
components, magnitudes, and phase angles) is included if the user requests the 
long form of the printed output. 
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6.2 Generalized Velocities 


In a manner similar to that for the generalized displacements of the proceeding 
section, the generalized velocities can be expressed in terms of the cos and sin 
components as follows: 

ZV^ = ZV^ cos (at + ZV* sin<at 
where: 

= generalized velocity for mode k 

By differentiating equation 67, we can express the amplitudes of the cos and sin 
components of the generalized velocities in terms of the generalized 
displacements found via the method of the proceeding section as follows: 


zn = < 

ZV- = -.az; 
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7.0 Equations For the Physical Quantities 

Once the generalized displacements and generalized velocities have been found, 
as outlined in section 6, they are used to calculate the physical quantities. 

The physical quantities are the last things calculated for each forced frequency 
solution. The following sections detail how the physical quantities are 
calculated. 

7.1 Physical Displacements, Velocities, and Modal Forces at 
the Points 


The physical displacement can be written in terms of the cos and sin components 
as follows: 

X.. = Xf . cos (ot + X*. sin u)t 
V tj y 


where: 


^ij = physical displacement for point i in direction j 

The amplitude of the cos and sin components of the physical displacement is 
found by summing over the modes as follows: 

n 

x^.= y <!>* zt 

y Ilk k 




« 

k 


where: 

^jk = displacement mode shape for point i, direction j, and mode k 

= amplitude of the cos component of the generalized displacement for 
mode k (see section 6.1) 

= amplitude of the sin component of the generalized displacement for 
mode k (see section 6.1) 

Similarly, the physical velocity can be written in terms of the cos and sin 
components as follows: 

V. . = V^. cos u>t + V*. sin bit 
y y y 

where: 

V.j = physical velocity for point i in direction j 
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The amplitude of the cos and sin components of the physical velocity is found by 
summing over the modes as follows: 

n 

V*. = y 

ij ^ ijk k 

k=l 


n 

V*. = y <D* zv! 

y ^ tjk k 
*=1 

where; 

ZV\ - amplitude of the cos component of the generalized velocity for mode k 
(see section 6.2) 

ZV* = amplitude of the sin component of the generalized velocity for mode k 
(see section 6.2) 

Similarly, the modal force can be written in terms of the cos and sin components 
as follows: 

F . . = F^.. cos (ot + F* . sin cot 
V y y 

where: 


^ij = modal force for point i in direction j 


The amplitude of the cos and sin components of the modal force is found by 
summing over the modes as follows: 


y 



c 

* 


F*. 

y 



where: 


= force mode shape for point i, direction j, and mode k 

Note that the physical displacements and physical velocities are calculated 
using the displacement mode shapes, while the modal forces are calculated using 
the force mode shapes. For the flexible vertical and horizontal plane 
subsystems, the displacement mode shapes are the translation and slope, while 
the force mode shapes are the shear and moment as entered on input sheet C-3. 
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So far we have expressed the physical displacements, physical velocities, and 
modal forces in terms of the cos and sin components. These quantities can also 
be expressed in terms of the magnitude and phase angle as follows: 

where: 



Usually, the user is primarily interested in the magnitudes of the quantities. 
For this reason, it is the magnitudes and phase angles of these quantities, 
rather than the cos and sin components, that are printed out and written to the 
output plot file. An exception to this is that the cos and sin components are 
also printed out if the long form of the printed output is requested via input 
variable lOUT on type A input sheet. 
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7.2 Physical Connecting and Gyroscopic Element Forces 


TETRA 2 has capability for six different types of physical connecting elements: 
type 1 (general spring-damper elements), type 2 (link elements), type 3 (rub 
elements), type 4 (engine support-links elements), type 5 (uncoupled point 
spring-damper elements), and type 6 (squeeze film damper elements). The first 
five types can be used for either transient or steady state analyses, while the 
type 6 (squeeze film damper elements) can be used only for transient analyses. 

In addition, TETRA 2 accounts for gyroscopic forces acting on a rotor by means 
of a gyroscopic element, which can be used for either transient or steady state 
analyses. Although not a physical connecting element, the gyroscopic element is 
classified as an element because gyroscopic forces are calculated very similarly 
to the damping forces of the physical connecting elements. This section 
concerns itself only with the element forces for steady state analysis runs, 
since the element forces for transient analysis runs were covered in reference 1 
and reference 3. 

The physical connecting or gyroscopic element force can be written in terms of 
the cos and sin components as follows: 

^ijk = ^jk 

where: 

^ijk = force that physical connecting or gyroscopic element k exerts on 
the engine components or ground for point i and direction j 

= amplitude of the cos component of the force that physical 
connecting or gyroscopic element k exerts on the engine component 
or ground for point i and direction j 

= amplitude of the sin component of the force that physical 
connecting or gyroscopic element k exerts on the engine component 
or ground for point i and direction j 

The amplitude of the cos and sin components of the force that the element exerts 
on the engine component or ground is calculated very differently for the 
nonlinear type 3 physical connecting element (rub element) than for the other 
elements. For the rub elements, these are calculated using harmonic averaging 
(except when the dead band has not been exceeded so that the rub elejnent forces 
are 0) as detailed in section 5. For joint I of the rub element, in the 

above equation is the same as the variables (for the vertical direction) 
and Ifor the horizontal direction) from section 5. Likewise, is the 

same as (for the vertical direction) and Ptm (for the horizontal 
direction) from section 5. The forces at joint J of the rub element are simply 
the negative of the forces at joint I of the rub element. 

For the other elements, on the other hand, the ijk and are calculated 
using the physical displacements and/or physical velocities at the joints of the 

element (the , x* , , and V*. that were found in section 7.1) and data 

y u y y 
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pertaining to the stiffness and damping of the element. Just what stiffness and 
damping data is used along with the physical displacements and/or physical 
velocities to calculate the element forces depends on the type of physical 
connecting element and the user input options which were chosen. Input 
stiffness matrix definition or input stiffness coefficients are used obtain the 
stiffness of the element. The damping of the element may be obtained by input 
damping matrix definition, input damping coefficients, or may be calculated 
using an input Q-factor and a frequency. The frequency used along with the 
input Q-factor to calculate the damping may either be input (non-structural 
damping), or the steady state forcing frequency or independent rotor speed may 
be used for this frequency (structural damping). In the case of the gyroscopic 
element, the polar moment of inertia and the rotor speed are used to calculate 
the damping. See reference 1 for more details about the physical connecting and 
gyroscopic elements. 

The physical connecting and gyroscopic element forces may also be expressed in 
terms of magnitude and phase angle as follows: 



Usually, the user is interested primarily in the magnitude of the element 
forces. For this reason, it is the magnitudes and phase angles of these 
quantities, rather than the cos and sin components, that are printed out and 
written to the output plot file. An exception to this is that the cos and sin 
components are also printed out if the long form of the printed output is 
requested via input variable lOUT on the type A input sheet. 
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7.3 Flexible Bladed Disk Displacements and Stresses 


For each flexible bladed disk, two modes are considered. These modes are the 
horizontal nodal diameter mode (referred to as mode P) and the vertical nodal 
diameter mode (referred to as mode Q). In a given TETRA model, there can be a 
maximum of two flexible bladed disks, and these flexible bladed disks must be 
located on the same rotor. See reference 2 for a detailed discussion of the two 
nodal diameter modes and flexible bladed disks in general. 


The equations for the displacements and stresses at a local point on flexible 
bladed disk number 1 or 2 (from reference 2 page 11 and 21) are: 


where: 


U= U (^P sinW + Q cosW^ 
V = V (^P sin 'P + Q cos 
S^= (^P sin 'P + Q cos 

Sj = S2 (^P sin V + Q cos 

S3 = (p sin 'P + Q cos 


( 68 ) 

(69) 

(70) 

(71) 

(72) 


U = tangential displacement of the local point on the flexible bladed disk 

U = input static (zero speed) mode shape for tangential translation of the 
local point on the flexible bladed disk 

V = axial displacement of the local point on the flexible bladed disk 


V = input static (zero speed) mode shape for axial translation of the 
local point on the flexible bladed disk 

= first stress component of the local point on the flexible bladed disk 

S^ = input modal stress for the first stress component of the local point 
on the flexible bladed disk 

= second stress component of the local point on the flexible bladed 
disk 

Sj • = input modal stress for the second stress component of the local point 
on the flexible bladed disk 

= third stress component of the local point on the flexible bladed disk 

S3 = input modal stress for the third stress component of the local point 
on the flexible bladed disk 
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P = generalized displacement for the horizontal nodal diameter mode of the 
flexible bladed disk 

Q = generalized displacement for the vertical nodal diameter mode of the 
flexible bladed disk 

V = polar angle of the local point on the flexible bladed 
disk 

The polar angle of the local point on the flexible bladed disk is 
found from: 

V= Qi+ <1> 


where; 

Q = flexible bladed disk rotor speed 
t = time 

<5 = input polar angle of the local point on the flexible bladed disk 
relative to the flexible bladed disk reference diameter (see input 
sheet C-15) 


Also, for a steady state analysis run we can express the generalized 
displacement of the P and Q modes in terms of magnitude and phase angle as 
follows: 

= P"* cos (cof - ^ J 
Q = Q"* cos (cot - 

where the magnitude and phase angle for the P and Q modes are calculated just 
like those of the other generalized displacements (see section 6.1). 


Plugging the expressions for 'P, P, and Q into equations 68 through 72, we arrive 
at fairly complex expressions for the displacements and stresses at a local 
point on a flexible bladed disk: 


U= U 


V= V 


Si=Si 


^2 


^ 3=53 


P"‘cos^ci>t — + Q"* cos ^Cl)^ — ^^^cos^Qf + <1* 

P"'cos^(ot — ^p^sin^Qt + + Q"* cos + <I> 

P"'cos^(ot - (ip^sin^Qt + <&^ + Q"* cos - <^jcos|^Q^ + <I> 

P”'cos^(ot - ^p^sin(^Qt + <I>j + Q"* cos - C,^coj(^ Qt + 4> 

P”‘cos^(o< — + Q*” cos — ^^cos(^Qt + «I> 
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We can define the magnitude portions of these expressions as 
follows: 


U"^ = 

p 

U"* = 

Q 


ytn 

p 

V"* = V Q 

9 


fj pm 

17 Q” 

y pm 
m 


cm . 

Qftl 


Sj P"* 

= s[q'^ 

S^pm 

S^Q"* 
= Fp'" 


^3,- 


S3Q" 


It is these magnitudes that are printed out and written onto the plot file for 
steady state analysis run. 

Substituting the magnitude quantities into the expressions for the displacement 
and stresses at a local point on a flexible bladed disk we get: 



8.0 CONCLUDING REMARKS FOR VOLUME 1. 


This volume documents the methodology for the steady state solution 
incorporated in TETRA 2. It is written to permit a straightforward 
understanding of its developments from the essentials of the theory, to the 
actual applications to engine dynamics and to the programmed working 
equations. These also include the treatment of nonlinear elements and 
the case of the intermittent rubs. 

It is intended that this volume should be a self contained description 
of the entire theory, as well as an accompaniment to the second volume. 

Volume 2 is the user's manual which contains both program input/output 
description and the trial or sample illustrative cases. The latter is a 
documentation of the progressive steps that were taken to debug and check 
the program, from simple degenerate cases to the twin spool engine model. 
This volume is also intended to be a self contained user's manual. However, 
to those interested in cross-checking program with theory. Volume 1 will 
be necessary. 
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